rm(list=ls())

# Package and data setup -----------------------------------------------------------

# Packages needed for this script
CRAN_packages <-
  c("tidyverse", "estimatr", "wesanderson", "colortools", "grid", "gtable", "stargazer",
    "patchwork", "ggstance", "haven", "cjoint")

# Install / Update packages available on CRAN 

# Packages already installed
installed_packages <- installed.packages()[,1]

# Check if you have the packages
check_packages <- CRAN_packages %in% installed_packages

# Install packages you don't have
not_installed <- CRAN_packages[!check_packages]
if (length(not_installed > 0)) {
  sapply(not_installed,install.packages)
}

# Load packages
library(tidyverse)
library(estimatr)
library(wesanderson)
library(colortools)
library(grid)
library(gtable)
library(stargazer)
library(patchwork)
library(ggstance)
library(haven)
library(cjoint)

# Remove unncessary objects
rm(check_packages, CRAN_packages, installed_packages, not_installed)


# Load data 
load("conjoint_data.RData")



# Figure 1 - How firms in sample resolved contract disputes ---------------
dispute_figures <- unique(conjoint_data %>% subset(dispute==1) %>% select(dispute, contains("dispute_res"), formal_firm, KEY)) 
dispute_figures$formal_or_informal <- ifelse(dispute_figures$formal_firm==1, "formal", "informal")

n_form <- length(dispute_figures$formal_or_informal[dispute_figures$formal_or_informal=="formal"])
n_inform <- length(dispute_figures$formal_or_informal[dispute_figures$formal_or_informal=="informal"])

amic_form <- table(dispute_figures$dispute_resolved_1, dispute_figures$formal_or_informal)[2,1]
amic_inform <- table(dispute_figures$dispute_resolved_1, dispute_figures$formal_or_informal)[2,2]

network_form <- table(dispute_figures$dispute_resolved_2, dispute_figures$formal_or_informal)[2,1]
network_inform <- table(dispute_figures$dispute_resolved_2, dispute_figures$formal_or_informal)[2,2]

gov_form <- table(dispute_figures$dispute_resolved_3, dispute_figures$formal_or_informal)[2,1]
gov_inform <- table(dispute_figures$dispute_resolved_3, dispute_figures$formal_or_informal)[2,2]

assoc_form <- table(dispute_figures$dispute_resolved_4, dispute_figures$formal_or_informal)[2,1]
assoc_inform <- table(dispute_figures$dispute_resolved_4, dispute_figures$formal_or_informal)[2,2]

rel_form <- table(dispute_figures$dispute_resolved_5, dispute_figures$formal_or_informal)[2,1]
rel_inform <- table(dispute_figures$dispute_resolved_5, dispute_figures$formal_or_informal)[2,2]

police_form <- table(dispute_figures$dispute_resolved_6, dispute_figures$formal_or_informal)[2,1]
police_inform <- table(dispute_figures$dispute_resolved_6, dispute_figures$formal_or_informal)[2,2]

lawyer_form <- table(dispute_figures$dispute_resolved_7, dispute_figures$formal_or_informal)[2,1]
lawyer_inform <- table(dispute_figures$dispute_resolved_7, dispute_figures$formal_or_informal)[2,2]

court_form <- table(dispute_figures$dispute_resolved_8, dispute_figures$formal_or_informal)[2,1]
court_inform <- table(dispute_figures$dispute_resolved_8, dispute_figures$formal_or_informal)[2,2]

never_form <- table(dispute_figures$dispute_resolved_9, dispute_figures$formal_or_informal)[2,1]
never_inform <- table(dispute_figures$dispute_resolved_9, dispute_figures$formal_or_informal)[2,2]

other_form <- table(dispute_figures$dispute_resolved_95, dispute_figures$formal_or_informal)[2,1]
other_inform <- table(dispute_figures$dispute_resolved_95, dispute_figures$formal_or_informal)[2,2]

dispute_type_n <- rep(c(1:9, 95),2)
dispute_numbers <- c(amic_form, network_form, gov_form, assoc_form, rel_form, police_form, lawyer_form, court_form , never_form, other_form , 
                     amic_inform, network_inform, gov_inform, assoc_inform, rel_inform, police_inform, lawyer_inform, court_inform , never_inform, other_inform )
dispute_pct <- c(((dispute_numbers[1:10]/n_form)*100), (dispute_numbers[11:20]/n_inform)*100)
firm_type <- c(rep("Formal firms", 10), rep("Informal firms", 10))
dispute_data_n <- data.frame(dispute_type_n, dispute_pct, firm_type)
dispute_data_n$dispute_type_n <- as.factor(dispute_data_n$dispute_type_n)
dispute_data_n$firm_type <- factor(dispute_data_n$firm_type, levels=c("Informal firms",
                                                                      "Formal firms"))
# Graph
dispute_data_n %>% ggplot(aes(dispute_type_n, dispute_pct)) +
  geom_bar(stat="identity") +
  facet_grid(~firm_type) +
  theme_bw() +
  scale_x_discrete(labels=c("Amicably",
                            "Social network",
                            "Gov. contact",
                            "Business assoc.",
                            "Religious group",
                            "Police",
                            "Lawyer",
                            "Court",
                            "Never resolved",
                            "Other")) +
  theme(axis.text = element_text( size = 12 ),
        axis.text.x = element_text( size = 14, angle = 45, hjust = 1 ),
        axis.title = element_text( size = 14, face = "bold" ),
        legend.position="none",
        strip.text = element_text(size = 14),
        plot.title = element_text(size = 16, face = "bold",hjust = 0.5)) +
  ylim(0,85) +
  ylab('Percentage by firm status') + xlab(label ="Dispute resolution method") + ggtitle('Dispute resolution by firm status')


# Figure 2: Distributions of firm size and valuation by formality  --------
fig2_data <- conjoint_data %>% select(firm_valuation, formal_firm, number_employees, KEY)
fig2_data <- unique(fig2_data)
fig2_data <- fig2_data %>% filter(firm_valuation!=98 & firm_valuation!=99)
fig2_data$formal_firm <- as.factor(fig2_data$formal_firm)
levels(fig2_data$formal_firm) <- c("Informal firms", "Formal firms")

# Firm valuation by formality status
fig2_data %>% ggplot(aes(as.factor(firm_valuation), group=formal_firm)) +
  geom_bar(aes(y = ..prop..), stat='count') + 
  scale_y_continuous(labels=scales::percent) +
  scale_x_discrete(labels=c("Less than \n 100.000",
                            "100.000- \n 500.000",
                            "500.001- \n 1.000.000",
                            "1.000.000- \n 5.000.000",
                            "5.000.001- \n 10.000.000",
                            "10.000.001- \n 100.000.000",
                            "More than \n 100.000.000")) +
  facet_grid(~formal_firm) + 
  theme_bw() + 
  theme(axis.text = element_text( size = 12 ),
        axis.text.x = element_text( size = 14, angle = 45, hjust = 1 ),
        axis.title = element_text( size = 14, face = "bold" ),
        legend.position="none",
        strip.text = element_text(size = 14),
        plot.title = element_text(size = 16, face = "bold",hjust = 0.5)) +
  ylab('Percentage by firm status') + xlab(label ="Firm valuation (FCFA)") + ggtitle('Firm valuation by formality status')

#Firm size by formality status
fig2_data %>% ggplot(aes(as.factor(number_employees), group=formal_firm)) +
  geom_bar(aes(y = ..prop..), stat='count') + 
  scale_y_continuous(labels=scales::percent) +
  scale_x_discrete(labels=c("1",
                            "2 to 5",
                            "6 to 10",
                            "11 to 20",
                            "21 to 50",
                            "More \n than 50")) +
  facet_grid(~formal_firm) + 
  theme_bw() + 
  theme(axis.text = element_text( size = 12 ),
        axis.text.x = element_text( size = 14 ),
        axis.title = element_text( size = 14, face = "bold" ),
        legend.position="none",
        strip.text = element_text(size = 14),
        plot.title = element_text(size = 16, face = "bold",hjust = 0.5)) +
  ylab('Percentage by firm status') + xlab(label ="Number of employees") + ggtitle('Firm size by formality status')


# Figure 3 - Main result: Influences of firms' likelihood to exchange --------------

# Regression output
reg1 <- lm_robust(formula = preferred_firm ~ c_contract + 
                    c_personal_cxn + c_political_affil + 
                    c_business_size + c_ethnicity + c_religion,
                  data = conjoint_data,
                  subset = conjoint_data$preferred_firm!=99 & conjoint_data$preferred_firm!=98,
                  clusters = KEY)
sum_reg1 <- summary(reg1)$coefficients

# Making data frame for figure
reg1_data_for_graph <- data.frame(Name = c(levels(conjoint_data$c_contract),
                                           levels(conjoint_data$c_personal_cxn), 
                                           levels(conjoint_data$c_political_affil),
                                           levels(conjoint_data$c_business_size),
                                           levels(conjoint_data$c_ethnicity),
                                           levels(conjoint_data$c_religion)
),
Factor = c(rep("Contract", 2), rep("Personal connections", 4), 
           rep("Political affiliation", 3), rep("Business size", 3),
           rep("Ethnicity", 5), rep("Religion", 5)
),
est = c(0, sum_reg1[2, 1], 0, sum_reg1[3:5, 1], 
        0, sum_reg1[6:7, 1], 0, sum_reg1[8:9, 1],
        0, sum_reg1[10:13, 1], 0, sum_reg1[14:17, 1]
),
se = c(0, sum_reg1[2, 2], 0, sum_reg1[3:5, 2], 
       0, sum_reg1[6:7, 2], 0, sum_reg1[8:9, 2],
       0, sum_reg1[10:13, 2], 0, sum_reg1[14:17, 2]
)
) %>%
  mutate(name = factor(Name, levels = rev( c(levels(conjoint_data$c_contract),
                                             levels(conjoint_data$c_personal_cxn), 
                                             levels(conjoint_data$c_political_affil),
                                             levels(conjoint_data$c_business_size),
                                             levels(conjoint_data$c_ethnicity),
                                             levels(conjoint_data$c_religion)
  )
  )
  ),
  attributes_ordered = factor(Factor, levels = unique(Factor)[c(5, 6, 3, 2, 4, 1)])
  )


# Figure
title <- expression(paste(Delta, " in Prob. of Respondent Choosing Firm, Relative to Reference Attribute"))
mtitle <- "AMCEs for Full Sample: Exchange Outcome"

ggplot(reg1_data_for_graph, aes(y = name, x = est, color = Factor)) +
  geom_vline(xintercept = 0, linetype = 2, color = "grey") + 
  facet_grid(attributes_ordered~., scale="free_y", space = "free_y") +
  geom_point(size = 2) + 
  geom_errorbarh(aes(xmin = est - 1.64 * se, xmax = est + 1.64 * se), height = 0, lwd = 1) +
  geom_errorbarh(aes(xmin = est - 1.96 * se, xmax = est + 1.96 * se), height = 0) +
  theme_bw() + theme(legend.position = "none") + 
  ylab(" ") + xlab(title) + ggtitle(mtitle) +
  scale_color_brewer(palette = "Dark2") +
  theme(axis.text.x = element_text(size = 14),
        axis.text.y = element_text(size = 14),
        strip.text = element_text(face="bold", size=10),
        text = element_text(size=14),
        legend.text=element_text(size=14))


# Figure 4 - Results conditional on co-ethnicity and co-religiosity -------

## 4a. Ethnicity ----

# Regressions
reg_coethnic <- lm_robust(formula = preferred_firm ~ c_contract + 
                            c_personal_cxn + c_political_affil + 
                            c_business_size + c_ethnicity + c_religion,
                          data = conjoint_data,
                          subset = conjoint_data$preferred_firm!=99 & conjoint_data$preferred_firm!=98 & conjoint_data$coethnic==1,
                          clusters = KEY)
sum_reg_coethnic <- summary(reg_coethnic)$coefficients

reg_noncoethnic <- lm_robust(formula = preferred_firm ~ c_contract + 
                               c_personal_cxn + c_political_affil + 
                               c_business_size + c_ethnicity + c_religion,
                             data = conjoint_data,
                             subset = conjoint_data$preferred_firm!=99 & conjoint_data$preferred_firm!=98 & conjoint_data$coethnic==0,
                             clusters = KEY)
sum_reg_noncoethnic <- summary(reg_noncoethnic)$coefficients

# Preparing graph data

#coethnic
reg_coethnic_data_for_graph <- data.frame(Name = c(levels(conjoint_data$c_ethnicity)),
                                          Factor = c(rep("Ethnicity", 5)),
                                          est = c(0,sum_reg_coethnic[10:13, 1]),
                                          se = c(0, sum_reg_coethnic[10:13, 2])) %>%
  mutate(name = factor(Name, levels = rev( c(levels(conjoint_data$c_ethnicity)))))

#noncoethnic
reg_noncoethnic_data_for_graph <- data.frame(Name = c(levels(conjoint_data$c_ethnicity)),
                                             Factor = c(rep("Ethnicity", 5)),
                                             est = c(0, sum_reg_noncoethnic[10:13, 1]),
                                             se = c(0, sum_reg_noncoethnic[10:13, 2])
) %>%
  mutate(name = factor(Name, levels = rev( c(levels(conjoint_data$c_ethnicity)))))

#combining above
reg_coethnic_data_for_graph$group <- "Coethnic"
reg_noncoethnic_data_for_graph$group <- "Non-coethnic"

combined_data_for_graph_ethnicity <- bind_rows(reg_coethnic_data_for_graph, reg_noncoethnic_data_for_graph)
combined_data_for_graph_ethnicity$group <- factor(combined_data_for_graph_ethnicity$group, levels = c("Coethnic", "Non-coethnic"))
combined_data_for_graph_ethnicity$name_order <- factor(combined_data_for_graph_ethnicity$name,
                                                       levels = rev(levels(combined_data_for_graph_ethnicity$name)))

legend_name <- ""

# Graph
coethnicity_graph <- ggplot(combined_data_for_graph_ethnicity, aes(y = name_order, x = est, color = group)) +
  geom_vline(xintercept = 0, linetype = 2, color = "grey") + 
  facet_grid(.~Factor, scale="free_y", space = "free_y") +
  geom_point(size = 2, position=ggstance::position_dodgev(height=0.2)) +
  geom_errorbarh(aes(xmin = est - 1.64 * se, xmax = est + 1.64 * se), 
                 height = 0, lwd = 1, position=ggstance::position_dodgev(height=0.2)) + 
  geom_errorbarh(aes(xmin = est - 1.96 * se, xmax = est + 1.96 * se), 
                 height = 0, position=ggstance::position_dodgev(height=0.2)) + 
  theme_bw() + 
  xlim(-.2, .2) +
  coord_flip() +
  theme(legend.position = "bottom") + 
  ylab(" ") + xlab(" ") + 
  scale_color_manual(legend_name, values=wes_palette(n=2, name="Royal1")) + 
  theme(axis.text.x = element_text(size = 14),
        axis.text.y = element_text(size = 14),
        strip.text = element_text(face="bold", size=10),
        text = element_text(size=14),
        legend.text=element_text(size=14))


## 4b. Religion ----

# Regressions
reg_coreligious <- lm_robust(formula = preferred_firm ~ c_contract + 
                               c_personal_cxn + c_political_affil + 
                               c_business_size + c_ethnicity + c_religion,
                             data = conjoint_data,
                             subset = conjoint_data$preferred_firm!=99 & conjoint_data$preferred_firm!=98 & conjoint_data$coreligion==1,
                             clusters = KEY)
sum_reg_coreligious <- summary(reg_coreligious)$coefficients


reg_noncoreligious <- lm_robust(formula = preferred_firm ~ c_contract + 
                                  c_personal_cxn + c_political_affil + 
                                  c_business_size + c_ethnicity + c_religion,
                                data = conjoint_data,
                                subset = conjoint_data$preferred_firm!=99 & conjoint_data$preferred_firm!=98 & conjoint_data$coreligion==0,
                                clusters = KEY)
sum_reg_noncoreligious <- summary(reg_noncoreligious)$coefficients

# Preparing graph data

# coreligious
reg_coreligious_data_for_graph <- data.frame(Name = c(levels(conjoint_data$c_religion)),
                                             Factor = c(rep("Religion", 5)),
                                             est = c(0, sum_reg_coreligious[14:17, 1]),
                                             se = c(0, sum_reg_coreligious[14:17, 2])) %>%
  mutate(name = factor(Name, levels = rev( c(
    levels(conjoint_data$c_religion)))))

# noncoreligious
reg_noncoreligious_data_for_graph <- data.frame(Name = c(levels(conjoint_data$c_religion)),
                                                Factor = c(rep("Religion", 5)),
                                                est = c(0, sum_reg_noncoreligious[14:17, 1]),
                                                se = c(0, sum_reg_noncoreligious[14:17, 2])) %>%
  mutate(name = factor(Name, levels = rev( c(levels(conjoint_data$c_religion)))))

#combining above
reg_coreligious_data_for_graph$group <- "Coreligious"
reg_noncoreligious_data_for_graph$group <- "Non-coreligious"
combined_data_for_graph_religion <- bind_rows(reg_coreligious_data_for_graph, reg_noncoreligious_data_for_graph)
combined_data_for_graph_religion$group <- factor(combined_data_for_graph_religion$group, levels = c("Coreligious", "Non-coreligious"))

legend_name <- ""


combined_data_for_graph_religion$name_ordered <- as.character(combined_data_for_graph_religion$name)
combined_data_for_graph_religion$name_ordered[combined_data_for_graph_religion$name_ordered=="Muslim (no brotherhood)"] <- "Muslim (no\nbrotherhood)"
combined_data_for_graph_religion$name_ordered <- factor(combined_data_for_graph_religion$name_ordered,
                                                        levels = combined_data_for_graph_religion$name_ordered[c(1:5)])


# Graph
coreligiousity_graph <- ggplot(combined_data_for_graph_religion, aes(y = name_ordered, x = est, color = group)) +
  geom_vline(xintercept = 0, linetype = 2, color = "grey") + 
  facet_grid(.~Factor, scale="free_y", space = "free_y") +
  geom_point(size = 2, position=ggstance::position_dodgev(height=0.2)) +
  geom_errorbarh(aes(xmin = est - 1.64 * se, xmax = est + 1.64 * se), 
                 height = 0, lwd = 1, position=ggstance::position_dodgev(height=0.2)) + 
  geom_errorbarh(aes(xmin = est - 1.96 * se, xmax = est + 1.96 * se), 
                 height = 0, position=ggstance::position_dodgev(height=0.2)) + 
  theme_bw() + 
  xlim(-.2, .2) +
  coord_flip() +
  theme(legend.position = "bottom") + 
  ylab(" ") + xlab(" ") + 
  scale_color_manual(legend_name, values=wes_palette(n=2, name="Royal1")) + 
  theme(axis.text.x = element_text(size = 14),
        axis.text.y = element_text(size = 14),
        strip.text = element_text(face="bold", size=10),
        text = element_text(size=14),
        legend.text=element_text(size=14))

# Combining into one figure:
coethnicity_graph | coreligiousity_graph



# Figure 5: Perceptions of likelihood of contract breach ------------------

# Regression output - prefer
reg1 <- lm_robust(formula = preferred_firm ~ c_contract + 
                    c_personal_cxn + c_political_affil + 
                    c_business_size + c_ethnicity + c_religion,
                  data = conjoint_data,
                  subset = conjoint_data$preferred_firm!=99 & conjoint_data$preferred_firm!=98,
                  clusters = KEY)
sum_reg1 <- summary(reg1)$coefficients


# Data for graph
reg1_data_for_graph <- data.frame(Name = c(levels(conjoint_data$c_contract),
                                           levels(conjoint_data$c_personal_cxn), 
                                           levels(conjoint_data$c_political_affil),
                                           levels(conjoint_data$c_business_size),
                                           levels(conjoint_data$c_ethnicity),
                                           levels(conjoint_data$c_religion)
),
Factor = c(rep("Contract", 2), rep("Personal connections", 4), 
           rep("Political affiliation", 3), rep("Business size", 3),
           rep("Ethnicity", 5), rep("Religion", 5)
),
est = c(0, sum_reg1[2, 1], 0, sum_reg1[3:5, 1], 
        0, sum_reg1[6:7, 1], 0, sum_reg1[8:9, 1],
        0, sum_reg1[10:13, 1], 0, sum_reg1[14:17, 1]
),
se = c(0, sum_reg1[2, 2], 0, sum_reg1[3:5, 2], 
       0, sum_reg1[6:7, 2], 0, sum_reg1[8:9, 2],
       0, sum_reg1[10:13, 2], 0, sum_reg1[14:17, 2]
)
) %>%
  mutate(name = factor(Name, levels = rev( c(levels(conjoint_data$c_contract),
                                             levels(conjoint_data$c_personal_cxn), 
                                             levels(conjoint_data$c_political_affil),
                                             levels(conjoint_data$c_business_size),
                                             levels(conjoint_data$c_ethnicity),
                                             levels(conjoint_data$c_religion)
  )
  )
  ),
  attributes_ordered = factor(Factor, levels = unique(Factor)[c(5, 6, 3, 2, 4, 1)])
  )

# Regression output - breach
reg2 <- lm_robust(formula = breach_risk ~ c_contract + 
                    c_personal_cxn + c_political_affil + 
                    c_business_size + c_ethnicity + c_religion,
                  data = conjoint_data,
                  subset = conjoint_data$breach_risk!=99 & conjoint_data$breach_risk!=98,
                  clusters = KEY)
sum_reg2 <- summary(reg2)$coefficients


# Data for graph
reg2_data_for_graph <- data.frame(Name = c(levels(conjoint_data$c_contract),
                                           levels(conjoint_data$c_personal_cxn), 
                                           levels(conjoint_data$c_political_affil),
                                           levels(conjoint_data$c_business_size),
                                           levels(conjoint_data$c_ethnicity),
                                           levels(conjoint_data$c_religion)),
                                  Factor = c(rep("Contract", 2), rep("Personal connections", 4), 
                                             rep("Political affiliation", 3), rep("Business size", 3),
                                             rep("Ethnicity", 5), rep("Religion", 5)
                                  ),
                                  est = c(0, sum_reg2[2, 1], 0, sum_reg2[3:5, 1], 
                                          0, sum_reg2[6:7, 1], 0, sum_reg2[8:9, 1],
                                          0, sum_reg2[10:13, 1], 0, sum_reg2[14:17, 1]
                                  ),
                                  se = c(0, sum_reg2[2, 2], 0, sum_reg2[3:5, 2], 
                                         0, sum_reg2[6:7, 2], 0, sum_reg2[8:9, 2],
                                         0, sum_reg2[10:13, 2], 0, sum_reg2[14:17, 2]
                                  )) %>%
  mutate(name = factor(Name, levels = rev( c(levels(conjoint_data$c_contract),
                                             levels(conjoint_data$c_personal_cxn), 
                                             levels(conjoint_data$c_political_affil),
                                             levels(conjoint_data$c_business_size),
                                             levels(conjoint_data$c_ethnicity),
                                             levels(conjoint_data$c_religion)
  )
  )
  ),
  attributes_ordered = factor(Factor, levels = unique(Factor)[c(5, 6, 3, 2, 4, 1)])
  )


reg1_data_for_graph$group <- "Exchange outcome"
reg2_data_for_graph$group <- "Breach outcome"

combined_data_for_graph <- bind_rows(reg1_data_for_graph, reg2_data_for_graph)
combined_data_for_graph$group <- factor(combined_data_for_graph$group, levels = c("Exchange outcome", "Breach outcome"))

title <- expression(paste(Delta, " in Prob. of Respondent Choosing Firm, Relative to Reference Attribute"))
legend_name <- ""
combined_title <- "AMCEs for Full Sample: Exchange and Breach Outcomes"

# Graph

ggplot(combined_data_for_graph, aes(y = name, x = est, color = group)) +
  geom_vline(xintercept = 0, linetype = 2, color = "grey") + 
  facet_grid(attributes_ordered~., scale="free_y", space = "free_y") +
  geom_point(size = 2, position=ggstance::position_dodgev(height=0.4)) +
  geom_errorbarh(aes(xmin = est - 1.64 * se, xmax = est + 1.64 * se), 
                 height = 0, lwd = 1, position=ggstance::position_dodgev(height=0.4)) + 
  geom_errorbarh(aes(xmin = est - 1.96 * se, xmax = est + 1.96 * se), 
                 height = 0, position=ggstance::position_dodgev(height=0.4)) + 
  theme_bw() + 
  theme(legend.position = "bottom") + 
  ylab(" ") + xlab(title) + ggtitle(combined_title) + 
  scale_color_manual(legend_name, values=wes_palette(n=2, name="Royal1")) +
  theme(axis.text.x = element_text(size = 14),
        axis.text.y = element_text(size = 14),
        strip.text = element_text(face="bold", size=10),
        text = element_text(size=14),
        legend.text=element_text(size=14))


# Figure 6 - Results conditional on firm formality status -------------------

conjoint_data$formal_firm <- as.factor(conjoint_data$formal_firm)

# Regression output
reg3 <- lm_robust(formula = preferred_firm ~ formal_firm*(c_contract + c_personal_cxn + c_political_affil + 
                                                            c_business_size + c_ethnicity + c_religion),
                  data = conjoint_data,
                  subset = conjoint_data$preferred_firm!=99 & conjoint_data$preferred_firm!=98,
                  clusters = KEY)
vcovmat <- reg3$vcov # variance covariance matrix
sum_reg3 <- summary(reg3)$coefficients #output of interest

g1 <- "Informal firms" 
g2 <- "Formal firms"
legend_name <- "" 

# Attribute labels for facet labels
attribute <-  c(
  rep("Contract\n[Ref: Informal\ncontract]", 1),
  rep("Personal connections\n[Ref: Not politically connected]", 3),
  rep("Political affiliation\n[Ref: No party affiliation]", 2), 
  rep("Business size\n[Ref: Medium business]", 2) ,
  rep("Ethnicity\n[Ref: Wolof]", 4), 
  rep("Religion\n[Ref: Muslim (no brotherhood)]", 4)
)

# Attribute levels (not including reference attributes)
trait <- c(
  c("Formal contract"), 
  c("Friend of local mayor", "Friend of MP", "Friend of president"),
  c("Opposition party member", "Ruling party member"),
  c("Large business", "Small business"),
  c("Diola", "Mandinuge", "Peul", "Serer"),
  c("Christian", "Layenne", "Mouride", "Tidjani")
)

trait[c(1, 2, 3, 4, 5, 6, 7, 8)] <- c("Formal\ncontract", "Knows\nlocal mayor", "Knows\nMP", "Knows\npresident", 
                                      "Opposition\nparty", "Ruling\nparty", "Large\nbusiness", "Small\nbusiness") 
group1 <- sum_reg3[3:18, 1:2] # coefs and SEs for group 1 (informal firms)
group2 <- cbind(sum_reg3[3:18, 1] + sum_reg3[19:34,1], # coefs and SEs for group 2 (formal firms)
                sqrt(diag(vcovmat)[3:18] + diag(vcovmat)[19:34] + 
                       2 * diag(vcovmat[3:18, 19:34])))
dif <- sum_reg3[19:34, 1:2] # coefs and SEs for difference
ests <- rbind(group1, group2, dif) 


hte_df <- data.frame(vars = rownames(ests), 
                     traits = rep(trait, 3), 
                     attributes = rep(attribute, 3), 
                     coef = ests[,1], 
                     se = ests[,2], 
                     Subsample = rep(c(g1, g2, "(Difference)"), each = 16),
                     panel = rep(c("AMCE Estimates", "AMCE Estimates", "Differences"), 
                                 each = 16)) %>%
  mutate(attributes_ordered = factor(attributes, levels = unique(attribute)[c(5, 6, 3, 2, 4, 1)]), # reorder attribute groupings as they'll appear in figure
         traits = factor(traits, level = trait),
         Subsample = factor(Subsample, level = c("Informal firms", "Formal firms", "(Difference)"))) #factor levels for legend 

title <- expression(paste(Delta, " in Prob. of Firm Being Chosen, Relative to Reference Attribute"))
mtitle <- "Conditional ACMEs by Firm Formality Status"

ggplot(hte_df, aes(y = coef, x = traits, color = Subsample))+
  geom_hline(yintercept = 0, linetype = 2, color = "grey") + 
  facet_grid(panel ~ attributes_ordered, scale="free_x", space = "free_x")+
  geom_point(size = 2, position = position_dodge(width = 0.3)) +
  geom_errorbar(aes(ymin = (coef - 1.96 * se), ymax = (coef+ 1.96 * se)), width=0,
                position = position_dodge(width = 0.3)) +
  geom_errorbar(aes(ymin = (coef - 1.64 * se), ymax = (coef+ 1.64 * se)), lwd=1, width = 0,
                position = position_dodge(width = 0.3)) +
  theme_bw() + xlab(" ")  + ylab(title) +
  scale_color_manual(legend_name, values=c(wes_palette(n=2, name="Royal1"), "gray55")) + 
  theme(axis.text.x = element_text(size = 10),
        axis.text.y = element_text(size = 10),
        strip.text = element_text(face="bold", size=8),
        text = element_text(size=14),
        legend.title=element_text(size=12),
        legend.text=element_text(size=12),
        legend.position = "bottom")



# Appendix ----------------------------------------------------------------


# Figure A1: Confidence in courts of law in Africa ------------------------

# Load Afrobarometer data directly from  website (internet connection required)
ab <- read_spss("https://afrobarometer.org/sites/default/files/data/round-6/merged_r6_data_2016_36countries2.sav",col_select = c("COUNTRY", "Q52J"))

# If Afrobarometer data is saved locally, instead executive the following commands:
#ab <- read_spss("merged_r6_data_2016_36countries2.sav")
#ab <- ab %>% dplyr::select("COUNTRY", "Q52J")

ab$name <- haven::as_factor(ab$COUNTRY)

# Removing non-answers/don't knows/etc.
ab$Q52J[ab$Q52J==-1 | ab$Q52J==9 | ab$Q52J==99] <- NA

means_by_country <- ab%>% group_by(name) %>% summarise(mean = mean(Q52J, na.rm=T))

# Ordering by outcome and producing graph
means_by_country <- means_by_country[order(means_by_country$mean),]
means_by_country$name <- factor(means_by_country$name, levels=means_by_country$name)
means_by_country <- means_by_country %>% mutate(ToHighlight = ifelse(name == "Senegal", "yes", "no"))

means_by_country %>% ggplot(aes(x=factor(name), y=mean, fill=ToHighlight)) +
  geom_bar(stat="identity") + coord_flip() + theme_bw() +
  ggtitle("") +
  theme(plot.title = element_text(hjust = 0.5)) +
  theme(plot.title = element_text(size=10)) +
  theme(axis.text=element_text(size=14),
        axis.title=element_text(size=14),
        title=element_text(size=14, face="bold"),
        legend.text=element_text(size=14)) +
  ylab("Trust in courts of law") +
  xlab("") +
  scale_y_continuous(limits = c(0,3), expand = c(0, 0)) +
  scale_fill_manual(values = c( "yes"="tomato", "no"="gray" ), guide = "none")


# Table B1: Sample summary statistics -------------------------------------

# Gathering variables of interest and subsetting to unique respondents
conjoint_data$formal_firm <- as.integer(conjoint_data$formal_firm)
data <- conjoint_data %>% select(KEY, 
                                 gender,
                                 age,
                                 educ,
                                 formal_firm,
                                 number_employees,
                                 monthly_revenue,
                                 firm_valuation,
                                 business_assoc,
                                 meet_other_businesses,
                                 access_credit_goodbad,
                                 declare_revenue,
                                 negotiate_taxes,
                                 taxes_paid,
                                 confidence_courts,
                                 dispute,
                                 percentage_formal,
                                 contract_more,
                                 contract_prefer,
                                 worked_state,
                                 polcon_useful,
                                 polcon_break,
                                 member_party,
                                 civil_assoc,
                                 contact_pol)
data <- unique(data)

# Re-categorizing the education variable to allow comparability
data <- data %>% mutate(educ_cat = ifelse(educ==0 | educ==1, 0, 
                                          ifelse(educ==2, 1,
                                                 ifelse(educ>40 & educ < 50, 2,
                                                        ifelse(educ>50 & educ<95, 3, educ)))))
# Summary function
sum_func <- function(variable){
  var_name_eval <- paste0("data$",variable)
  var_name <- eval(expr = parse(text = var_name_eval))
  other <- length(var_name[var_name==95])
  refuse <- length(var_name[var_name==98])
  missing <- length(var_name[var_name==99])
  sum_eval <- paste0("data %>% select(",variable,") %>% filter(",variable,"<95) %>% summarise_all(funs(mean = mean, sd = sd, min = min, median = median, max = max))")
  sums <- eval(expr = parse(text = sum_eval))
  full_sum_eval <- paste0("cbind('",variable,"',sums, other, refuse, missing)")
  full_sum <- eval(expr = parse(text = full_sum_eval))
  names(full_sum)[1] <- "variable_name"
  return(full_sum)
} 



#variables to summarize
var.list <- c("gender",
              "age",
              "educ_cat",
              "formal_firm",
              "number_employees",
              "monthly_revenue",
              "firm_valuation",
              "business_assoc",
              "meet_other_businesses",
              "access_credit_goodbad",
              "declare_revenue",
              "negotiate_taxes",
              "taxes_paid",
              "confidence_courts",
              "dispute",
              "percentage_formal",
              "contract_more",
              "contract_prefer",
              "worked_state",
              "polcon_useful",
              "polcon_break",
              "member_party",
              "civil_assoc",
              "contact_pol")

# Applying function to variables of interest and generating summary table
summing <- as.data.frame(lapply(var.list, sum_func)) 
sum_tab <- as.data.frame(matrix(summing, ncol=9, nrow=(length(summing)/9),byrow = T)) 
sum_tab <- apply(sum_tab[2:9],MARGIN = 2,FUN =  as.numeric) 
sum_tab <- data.frame(var.list, sum_tab) 
names(sum_tab) <- c("Variable", "Mean", "Std. Dev.", "Min", "Median", "Max", "Other", "Refusal", "Don't know")
stargazer(sum_tab, summary=F, digits = 3, float=F, rownames = F)


# Appendix F: Diagnostic tests -----------------------

# Figure F3: No carryover effects 

conjoint_nonas_prefer <- conjoint_data %>% filter(prefer!=98 & prefer!=99)
conjoint_nonas_prefer$round <- as.factor(conjoint_nonas_prefer$round)


amce_prefer_carryover <- amce(formula = preferred_firm ~ 
                                `c_ethnicity`*round + 
                                `c_religion`*round +
                                `c_political_affil`*round + 
                                `c_personal_cxn`*round + 
                                `c_business_size`*round +
                                `c_contract`*round,
                              na.ignore=TRUE,
                              respondent.varying = "round",
                              cluster = TRUE,
                              respondent.id = "KEY",
                              data = conjoint_nonas_prefer)


attribute_names <- c("Business size", 
                     "Contract", 
                     "Ethnicity", 
                     "Personal connections",
                     "Political affiliation", 
                     "Religion")

plot(amce_prefer_carryover, main = "Outcome: Which deal are you likely to choose",
     attribute.names = attribute_names,
     plot.display="interaction")


# Figure F4: No profile-order effects

conjoint_nonas_prefer <- conjoint_data %>% filter(prefer!=98 & prefer!=99)
conjoint_nonas_prefer$profile <- as.factor(conjoint_nonas_prefer$profile)

amce_prefer_profile_order <- amce(formula = preferred_firm ~ 
                                    `c_ethnicity`*profile + 
                                    `c_religion`*profile +
                                    `c_political_affil`*profile + 
                                    `c_personal_cxn`*profile + 
                                    `c_business_size`*profile +
                                    `c_contract`*profile,
                                  na.ignore=TRUE,
                                  respondent.varying = "profile",
                                  cluster = TRUE,
                                  respondent.id = "KEY",
                                  data = conjoint_nonas_prefer)

attribute_names <- c("Business size", 
                     "Contract", 
                     "Ethnicity", 
                     "Personal connections",
                     "Political affiliation", 
                     "Religion")

plot(amce_prefer_profile_order, main = "Outcome: Which deal are you likely to choose",
     attribute.names = attribute_names,
     plot.display="interaction")


# Table F2: Balance across treatment groups 
round((table(conjoint_data$c_political_affil)/nrow(conjoint_data))*100, 1)
round((table(conjoint_data$c_business_size)/nrow(conjoint_data))*100, 1)
round((table(conjoint_data$c_contract)/nrow(conjoint_data))*100, 1)
round((table(conjoint_data$c_personal_cxn)/nrow(conjoint_data))*100, 1)
round((table(conjoint_data$c_ethnicity)/nrow(conjoint_data))*100, 1)
round((table(conjoint_data$c_religion)/nrow(conjoint_data))*100, 1)


# Table F3: Balance across respondent characteristics
vars <- conjoint_data %>% select(confidence_courts, worked_state, gender, formal_firm, number_employees,
                                 monthly_revenue, firm_valuation,business_formalized_after_start,
                                 business_assoc, dispute, 
                                 c_ethnicity, c_religion, c_political_affil, c_personal_cxn,
                                 c_business_size, c_contract, KEY)

#clustered std error function
vcovCluster <- function(model,cluster)
{
  require(sandwich)
  require(lmtest)
  if(nrow(model.matrix(model))!=length(cluster)){
    stop("check your data: cluster variable has different N than model")
  }
  M <- length(unique(cluster))
  N <- length(cluster)           
  K <- model$rank   
  if(M<50){
    warning("Fewer than 50 clusters, variances may be unreliable (could try block bootstrap instead).")
  }
  dfc <- (M/(M - 1)) * ((N - 1)/(N - K))
  uj  <- apply(estfun(model), 2, function(x) tapply(x, cluster, sum));
  rcse.cov <- dfc * sandwich(model, meat = crossprod(uj)/N)
  return(rcse.cov)
}

#function that outputs models with standard errors clustered by respondent
balance_conjoint <- function(x){
  ols_model <- lm(x ~ as.factor(c_ethnicity) +
                    as.factor(c_religion) +
                    as.factor(c_political_affil) +
                    as.factor(c_personal_cxn) +
                    as.factor(c_business_size) +
                    as.factor(c_contract), 
                  data = vars)
  coeftest(ols_model, vcovCluster(ols_model, cluster=vars$KEY))
}

#imputing means for don't know/refuse to respond responses
vars$confidence_courts[vars$confidence_courts>90] <- NA 
vars$confidence_courts[is.na(vars$confidence_courts)==T] <- mean(vars$confidence_courts, na.rm=T)

vars$monthly_revenue[vars$monthly_revenue>90] <- NA 
vars$monthly_revenue[is.na(vars$monthly_revenue)==T] <- mean(vars$monthly_revenue, na.rm=T)

vars$firm_valuation[vars$firm_valuation>90] <- NA 
vars$firm_valuation[is.na(vars$firm_valuation)==T] <- mean(vars$firm_valuation, na.rm=T)

vars$business_formalized_after_start[vars$business_formalized_after_start>90] <- NA 
vars$business_formalized_after_start[is.na(vars$business_formalized_after_start)==T] <- mean(vars$business_formalized_after_start, na.rm=T)

#balance table
stargazer(lapply(X = vars[1:10],FUN = balance_conjoint), float=F)


# Figure F5: No attribute order effects

conjoint_nonas_prefer <- conjoint_data %>% filter(prefer!=98 & prefer!=99)

att_order <- amce(formula = preferred_firm ~ 
                    `c_ethnicity`*`Randomization scheme` + 
                    `c_religion`*`Randomization scheme` +
                    `c_political_affil`*`Randomization scheme` + 
                    `c_personal_cxn`*`Randomization scheme` + 
                    `c_business_size`*`Randomization scheme` +
                    `c_contract`*`Randomization scheme`,
                  na.ignore=TRUE,
                  respondent.varying = "`Randomization scheme`",
                  cluster = TRUE,
                  respondent.id = "KEY",
                  data = conjoint_nonas_prefer)

attribute_names <- c("Business size", 
                     "Contract", 
                     "Ethnicity", 
                     "Personal connections",
                     "Political affiliation", 
                     "Religion")


plot(att_order, main = "Outcome: Which deal are you likely to choose",
     attribute.names = attribute_names,
     plot.display="interaction")

# Figure F6: Attrition as predicted by treatment
conjoint_data <- conjoint_data %>% mutate(attritted = 
                                            ifelse(prefer==98|prefer==99, 1, 0))

amce_attrit <- amce(formula = attritted ~ c_contract + 
                      c_personal_cxn + c_political_affil + 
                      c_business_size + c_ethnicity + c_religion,
                    na.ignore=TRUE,
                    cluster = TRUE,
                    respondent.id = "KEY",
                    data = conjoint_data)

attribute_names <- c("Business size", 
                     "Contract", 
                     "Ethnicity", 
                     "Personal connections",
                     "Political affiliation", 
                     "Religion")

plot(amce_attrit, main = "Outcome: Attrition",
     attribute.names = attribute_names)

# Appendix G: Corresponding tables for figure results ---------------------

# Table G5: AMCE results for prefer outcome
reg1 <- lm_robust(formula = preferred_firm ~ c_contract + 
                    c_personal_cxn + c_political_affil + 
                    c_business_size + c_ethnicity + c_religion,
                  data = conjoint_data,
                  subset = conjoint_data$preferred_firm!=99 & conjoint_data$preferred_firm!=98,
                  clusters = KEY)
summary(reg1) 

# Table G6: AMCE results for breach outcome
reg2 <- lm_robust(formula = breach_risk ~ c_contract + 
                    c_personal_cxn + c_political_affil + 
                    c_business_size + c_ethnicity + c_religion,
                  data = conjoint_data,
                  subset = conjoint_data$breach_risk!=99 & conjoint_data$breach_risk!=98,
                  clusters = KEY)
summary(reg2)


# Table G7: AMCE results by firm formality
reg3 <- lm_robust(formula = preferred_firm ~ formal_firm*(c_contract + c_personal_cxn + c_political_affil + 
                                                            c_business_size + c_ethnicity + c_religion),
                  data = conjoint_data,
                  subset = conjoint_data$preferred_firm!=99 & conjoint_data$preferred_firm!=98,
                  clusters = KEY)
summary(reg3)

# Table G8: AMCE results by co-ethnicity
reg_coethnic <- lm_robust(formula = preferred_firm ~ coethnic*(c_contract + 
                                                                 c_personal_cxn + c_political_affil + 
                                                                 c_business_size + c_ethnicity + c_religion),
                          data = conjoint_data,
                          subset = conjoint_data$preferred_firm!=99 & conjoint_data$preferred_firm!=98,
                          clusters = KEY)
summary(reg_coethnic)


# Table G9: AMCE results by co-religion
reg_coreligious <- lm_robust(formula = preferred_firm ~ coreligion*(c_contract + 
                                                                      c_personal_cxn + c_political_affil + 
                                                                      c_business_size + c_ethnicity + c_religion),
                             data = conjoint_data,
                             subset = conjoint_data$preferred_firm!=99 & conjoint_data$preferred_firm!=98,
                             clusters = KEY)
summary(reg_coreligious)

# Appendix H: Effects by respondents’ political affiliations --------------
conjoint_nonas_prefer <- conjoint_data %>% filter(prefer!=98 & prefer!=99)


conjoint_nonas_prefer <- conjoint_nonas_prefer %>% mutate(ruling_or_opp = 
                                                            ifelse(party==1 | party ==5 | 
                                                                     party_other=="Benno Bokk Yakaar" |
                                                                     party_other=="Apr", 
                                                                   "ruling", "opposition"))

conjoint_nonas_prefer <- conjoint_nonas_prefer %>% mutate(Party = case_when(
  member_party==0 ~ "No party",
  ruling_or_opp== "ruling" ~ "Ruling party respondents",
  ruling_or_opp== "opposition" ~ "Opposition party respondents",
  TRUE ~ "NA"
))

conjoint_nonas_prefer$Party[conjoint_nonas_prefer$Party=="NA"] <- "No party"

conjoint_nonas_prefer$Party <- as.factor(conjoint_nonas_prefer$Party)

distaste <- amce(formula = preferred_firm ~ c_contract + 
                   c_personal_cxn + c_political_affil*Party + 
                   c_business_size + c_ethnicity + c_religion,
                 na.ignore=TRUE,
                 cluster = TRUE,
                 respondent.varying = "Party",
                 respondent.id = "KEY",
                 data = conjoint_nonas_prefer)


plot(distaste, main = "Outcome: Which deal are you likely to choose",
     plot.display="interaction"
)

# Appendix I: ACIEs -------------------------------------------------------
reg_acie <- lm_robust(formula = preferred_firm ~ c_contract*c_political_affil + 
                        c_contract*c_personal_cxn + c_political_affil*c_personal_cxn  + 
                        c_contract*c_business_size + c_contract*c_ethnicity + c_contract*c_religion,
                      data = conjoint_data,
                      subset = conjoint_data$preferred_firm!=99 & conjoint_data$preferred_firm!=98,
                      clusters = KEY)
summary(reg_acie)



# Appendix J: Standard errors clustered by district -----------------------


# Regression output - prefer
reg1 <- lm_robust(formula = preferred_firm ~ c_contract + 
                    c_personal_cxn + c_political_affil + 
                    c_business_size + c_ethnicity + c_religion,
                  data = conjoint_data,
                  subset = conjoint_data$preferred_firm!=99 & conjoint_data$preferred_firm!=98,
                  clusters = neighborhood)
sum_reg1 <- summary(reg1)$coefficients


# Data for graph
reg1_data_for_graph <- data.frame(Name = c(levels(conjoint_data$c_contract),
                                           levels(conjoint_data$c_personal_cxn), 
                                           levels(conjoint_data$c_political_affil),
                                           levels(conjoint_data$c_business_size),
                                           levels(conjoint_data$c_ethnicity),
                                           levels(conjoint_data$c_religion)
),
Factor = c(rep("Contract", 2), rep("Personal connections", 4), 
           rep("Political affiliation", 3), rep("Business size", 3),
           rep("Ethnicity", 5), rep("Religion", 5)
),
est = c(0, sum_reg1[2, 1], 0, sum_reg1[3:5, 1], 
        0, sum_reg1[6:7, 1], 0, sum_reg1[8:9, 1],
        0, sum_reg1[10:13, 1], 0, sum_reg1[14:17, 1]
),
se = c(0, sum_reg1[2, 2], 0, sum_reg1[3:5, 2], 
       0, sum_reg1[6:7, 2], 0, sum_reg1[8:9, 2],
       0, sum_reg1[10:13, 2], 0, sum_reg1[14:17, 2]
)
) %>%
  mutate(name = factor(Name, levels = rev( c(levels(conjoint_data$c_contract),
                                             levels(conjoint_data$c_personal_cxn), 
                                             levels(conjoint_data$c_political_affil),
                                             levels(conjoint_data$c_business_size),
                                             levels(conjoint_data$c_ethnicity),
                                             levels(conjoint_data$c_religion)
  )
  )
  ),
  attributes_ordered = factor(Factor, levels = unique(Factor)[c(5, 6, 3, 2, 4, 1)])
  )

# Regression output - breach
reg2 <- lm_robust(formula = breach_risk ~ c_contract + 
                    c_personal_cxn + c_political_affil + 
                    c_business_size + c_ethnicity + c_religion,
                  data = conjoint_data,
                  subset = conjoint_data$breach_risk!=99 & conjoint_data$breach_risk!=98,
                  clusters = neighborhood)
sum_reg2 <- summary(reg2)$coefficients


# Data for graph
reg2_data_for_graph <- data.frame(Name = c(levels(conjoint_data$c_contract),
                                           levels(conjoint_data$c_personal_cxn), 
                                           levels(conjoint_data$c_political_affil),
                                           levels(conjoint_data$c_business_size),
                                           levels(conjoint_data$c_ethnicity),
                                           levels(conjoint_data$c_religion)),
                                  Factor = c(rep("Contract", 2), rep("Personal connections", 4), 
                                             rep("Political affiliation", 3), rep("Business size", 3),
                                             rep("Ethnicity", 5), rep("Religion", 5)
                                  ),
                                  est = c(0, sum_reg2[2, 1], 0, sum_reg2[3:5, 1], 
                                          0, sum_reg2[6:7, 1], 0, sum_reg2[8:9, 1],
                                          0, sum_reg2[10:13, 1], 0, sum_reg2[14:17, 1]
                                  ),
                                  se = c(0, sum_reg2[2, 2], 0, sum_reg2[3:5, 2], 
                                         0, sum_reg2[6:7, 2], 0, sum_reg2[8:9, 2],
                                         0, sum_reg2[10:13, 2], 0, sum_reg2[14:17, 2]
                                  )) %>%
  mutate(name = factor(Name, levels = rev( c(levels(conjoint_data$c_contract),
                                             levels(conjoint_data$c_personal_cxn), 
                                             levels(conjoint_data$c_political_affil),
                                             levels(conjoint_data$c_business_size),
                                             levels(conjoint_data$c_ethnicity),
                                             levels(conjoint_data$c_religion)
  )
  )
  ),
  attributes_ordered = factor(Factor, levels = unique(Factor)[c(5, 6, 3, 2, 4, 1)])
  )


reg1_data_for_graph$group <- "Exchange outcome"
reg2_data_for_graph$group <- "Breach outcome"

combined_data_for_graph <- bind_rows(reg1_data_for_graph, reg2_data_for_graph)
combined_data_for_graph$group <- factor(combined_data_for_graph$group, levels = c("Exchange outcome", "Breach outcome"))

title <- expression(paste(Delta, " in Prob. of Respondent Choosing Firm, Relative to Reference Attribute"))
legend_name <- ""
combined_title <- "AMCEs for Full Sample: Exchange and Breach Outcomes"

# Graph

ggplot(combined_data_for_graph, aes(y = name, x = est, color = group)) +
  geom_vline(xintercept = 0, linetype = 2, color = "grey") + 
  facet_grid(attributes_ordered~., scale="free_y", space = "free_y") +
  geom_point(size = 2, position=ggstance::position_dodgev(height=0.4)) +
  geom_errorbarh(aes(xmin = est - 1.64 * se, xmax = est + 1.64 * se), 
                 height = 0, lwd = 1, position=ggstance::position_dodgev(height=0.4)) + 
  geom_errorbarh(aes(xmin = est - 1.96 * se, xmax = est + 1.96 * se), 
                 height = 0, position=ggstance::position_dodgev(height=0.4)) + 
  theme_bw() + 
  theme(legend.position = "bottom") + 
  ylab(" ") + xlab(title) + ggtitle(combined_title) + 
  scale_color_manual(legend_name, values=wes_palette(n=2, name="Royal1")) +
  theme(axis.text.x = element_text(size = 14),
        axis.text.y = element_text(size = 14),
        strip.text = element_text(face="bold", size=10),
        text = element_text(size=14),
        legend.text=element_text(size=14))



